####### Michael Shepherd
###### Politics of Rural Hospital Closures 
######################################################################
#                   Loading Packages 
######################################################################
rm(list = ls()) 
library(tidyverse) 
library(ggplot2)
library(readr)
library(ggpubr)
library(sandwich)
library(lmtest)
library(stargazer)
library(haven)
library(readxl)
library(MatchIt)

setwd("~/Desktop/RURALCLOSURES")

######################################################################
#     Loading, Manipulating, and Merging County/ Hospital Data
######################################################################
# County RUCA Codes 
rural <- read_csv("Rural-urban Continuum Code 2013-Table 1.csv")

# Non-metro as rural 
rural$rural <- NA
rural$rural[rural$RUCC_2013<=3] <- 0
rural$rural[rural$RUCC_2013>=4] <- 1

# Loading Hospital Closure Data 
hosp <- read_csv("hospital_closure_counties_FINAL.csv")
# copying to perserve original 
hosp2 <- hosp 
hosp2$Closure <- 1 # Closure Indicator 

hosp3 <- hosp2 %>%
  group_by(FIPS,`Closure Year`) %>%
  summarise(Closure = sum(Closure, na.rm=T)) # Closures by County Year  
# Milan County, TX Lost Two Hospitals in 2018 (dropped in main analyses but reported in robusts). 

hosp4 <- merge(rural, hosp3, by="FIPS", all=T)

# Creating Treatment Indicator in Main Text 
hosp4$treatedever <- NA 
hosp4$treatedever[hosp4$Closure==1] <- 1 
hosp4$treatedever[is.na(hosp4$Closure)] <- 0

# Creating Treatment Indicator Including Texas County losing Two in One Year
hosp4$treatedever2 <- NA 
hosp4$treatedever2[hosp4$Closure>=1] <- 1 
hosp4$treatedever2[is.na(hosp4$Closure)] <- 0

# RWJ Health Data for Health Matching Controls 
rwj <- read_csv("rwj_2011_state_county.csv")
rwj$FIPS <- rwj$`5-digit FIPS Code`

#Location of Open Hospitals in US 
hospitals2018 <- read_excel("AllHospitalList2018.xlsx")

hospitals2018$FIPS <- as.numeric(hospitals2018$FIPS)

hospitals2018$Hospitals <- 1 

open <-hospitals2018 %>%
  group_by(FIPS) %>%
  summarise(Hospitals= sum(Hospitals, na.rm=T)) # Number of Hospitals in counties with hospitals 

open2 <-hospitals2018 %>%
  group_by(STATE, FIPS) %>%
  summarise(Hospitals= sum(Hospitals, na.rm=T))# Number of Hospitals in counties with hospitals by state


hosp4 <- merge(hosp4, rwj, by="FIPS")
hosp4$treatedever[is.na(hosp4$treatedever)] <- 0 # 0 For Counties Never Treated 
hosp4$treatedever2[is.na(hosp4$treatedever2)] <- 0 # 0 For Counties Never Treated (Robust to Milam)

# Subsetting to Analysis Set 
hosp4 <- hosp4[which(hosp4$State!="CA" | hosp4$State!="WA" |  hosp4$State!="RI" |  hosp4$State!="NY"| hosp4$State!="MN"),]
hosp4 <- hosp4[which(hosp4$`Closure Year`>=2008 | hosp4$treatedever==0),]
hosp6 <- merge(hosp4, open, by="FIPS", all=T)
hosp6$Hospitals[is.na(hosp6$Hospitals)] <- 0 # Missing equals no hospitals 
hosp6 <- hosp6[which(hosp6$State!="CA"),] # Be sure Blue State Closures/PotentialMatches removed before final matching
hosp6 <- hosp6[which(hosp6$State!="WA"),]
hosp6 <- hosp6[which(hosp6$State!="RI"),]
hosp6 <- hosp6[which(hosp6$State!="NY"),]
hosp6 <- hosp6[which(hosp6$State!="MN"),]

                     
hosp6r <- hosp6 # Milam TX Robustness Set
hosp6r$treatedever2[is.na(hosp6r$treatedever2)] <- 0  # Missing equals not treated

hosp6 <- hosp6[which(hosp6$Hospitals<=2),] # Counties without numerous hospitals for primary matching
hosp6$treatedever[is.na(hosp6$treatedever)] <- 0 # Missing equals not treated


# Subset to complete cases for Matching 
hosp6 <- hosp6[complete.cases(hosp6$`Uninsured adults raw value`), ] 
hosp6 <- hosp6[complete.cases(hosp6$`Median household income raw value`), ] 
hosp6 <- hosp6[complete.cases(hosp6$`Poor or fair health raw value`), ] 
hosp6 <- hosp6[complete.cases(hosp6$`% Rural raw value`), ] 


# Subset to complete cases for Matching plus Milam Robust
hosp6r <- hosp6r[complete.cases(hosp6r$`Uninsured adults raw value`), ] 
hosp6r <- hosp6r[complete.cases(hosp6r$`Median household income raw value`), ] 
hosp6r <- hosp6r[complete.cases(hosp6r$`Poor or fair health raw value`), ] 
hosp6r <- hosp6r[complete.cases(hosp6r$`% Rural raw value`), ] 



# Non-expansion status by 2018
hosp6$nonexpansion <- 0 
hosp6$nonexpansion[hosp6$State=="TX" |
                     hosp6$State=="TN" |
                     hosp6$State=="FL" |
                     hosp6$State=="NC" |
                     hosp6$State=="GA" |  hosp6$State=="AL" | hosp6$State=="MS" | 
                     hosp6$State=="KS" | hosp6$State=="WY" | hosp6$State=="WI" | 
                     hosp6$State=="SC" | hosp6$State=="SD" | hosp6$State=="OK" | hosp6$State=="MO" | hosp6$State=="NE" | hosp6$State=="ID" |
                     hosp6$State=="ME"] <- 1 


# Non-expansion status by 2018 (for Milam Robust set)
hosp6r$nonexpansion <- 0 
hosp6r$nonexpansion[hosp6r$State=="TX" |
                     hosp6r$State=="TN" |
                     hosp6r$State=="FL" |
                     hosp6r$State=="NC" |
                     hosp6r$State=="GA" |  hosp6r$State=="AL" | hosp6r$State=="MS" | 
                     hosp6r$State=="KS" | hosp6r$State=="WY" | hosp6r$State=="WI" | 
                     hosp6r$State=="SC" | hosp6r$State=="SD" | hosp6r$State=="OK" | hosp6r$State=="MO" | hosp6r$State=="NE" | hosp6r$State=="ID" |
                     hosp6r$State=="ME"] <- 1 



######################################################################
#               Matching Processes  
######################################################################

# Nearest Neighbor Matchingng Used in Main Analyses in Article (no Milam County)
m.out0 <- matchit(treatedever ~ `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
                    `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + nonexpansion, data = hosp6,
                  method = "nearest", distance = "glm")
summary(m.out0, un = FALSE)
plot(m.out0, type = "jitter", interactive = FALSE)
plot(summary(m.out0))

m.data <- match.data(m.out0)


# Nearest Neighbor Matching Primary Strategy Including Milam County 
m.out0r <- matchit(treatedever2 ~ `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
                    `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + nonexpansion, data = hosp6r,
                  method = "nearest", distance = "glm")
summary(m.out0r, un = FALSE)
plot(m.out0r, type = "jitter", interactive = FALSE)
plot(summary(m.out0r))


m.datar <- match.data(m.out0r)


# Nearest Neighbor Matching Exact Within State Matching Robustness Check 
m.out0_state <- matchit(treatedever ~ `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
                          `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + nonexpansion, data = hosp6,
                        method = "nearest", exact = "State", distance = "glm")
summary(m.out0_state, un = FALSE)
plot(m.out0_state, type = "jitter", interactive = FALSE)
plot(summary(m.out0_state))

m.out0_statematch <- match.data(m.out0_state)



# Nearest Neighbor Matching Exact Within State Matching Robustness Check with Milam County, TX 
m.out0rs <- matchit(treatedever2 ~ `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
                     `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + nonexpansion, data = hosp6r,
                   method = "nearest", , exact = "State", distance = "glm")
summary(m.out0rs, un = FALSE)
plot(m.out0rs, type = "jitter", interactive = FALSE)
plot(summary(m.out0rs))


m.datars <- match.data(m.out0rs)


######################################################################
#           Probing Geographic Balance for All Matched Sets
######################################################################


matchedhospitalcounties <- m.data
matchedhospitalcountiesr <- m.datar
matchedhospitalcountiesrns <- m.datars
matchedhospitalcountiesns <- m.out0_statematch

# Observational Balance Post-Matching (Primary Matched Data)
summary(lm(matchedhospitalcounties$`% Rural raw value` ~ matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$`% Non-Hispanic African American raw value` ~ matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$`% Hispanic raw value` ~ matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$`% 65 and older raw value` ~ matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$`Median household income raw value` ~ matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$`Poor or fair health raw value`~ matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$`Uninsured adults raw value`~ matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$Hospitals~ matchedhospitalcounties$treatedever))

# Observational Balance Post-Matching (Primary Matched + Texas Case )
summary(lm(matchedhospitalcountiesr$`% Rural raw value` ~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$`% Non-Hispanic African American raw value` ~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$`% Hispanic raw value` ~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$`% 65 and older raw value` ~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$`Median household income raw value` ~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$`Poor or fair health raw value`~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$`Uninsured adults raw value`~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$Hospitals~ matchedhospitalcountiesr$treatedever2))



# Observational Balance Post-Matching (In State Matched )
summary(lm(matchedhospitalcountiesns$`% Rural raw value` ~ matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$`% Non-Hispanic African American raw value` ~ matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$`% Hispanic raw value` ~ matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$`% 65 and older raw value` ~matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$`Median household income raw value` ~ matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$`Poor or fair health raw value`~ matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$`Uninsured adults raw value`~ matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$Hospitals~ matchedhospitalcountiesns$treatedever))

# Observational Balance Post-Matching (In State Matched + Texas Case )
summary(lm(matchedhospitalcountiesrns$`% Rural raw value` ~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$`% Non-Hispanic African American raw value` ~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$`% Hispanic raw value` ~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$`% 65 and older raw value` ~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$`Median household income raw value` ~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$`Poor or fair health raw value`~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$`Uninsured adults raw value`~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$Hospitals~ matchedhospitalcountiesrns$treatedever2))


# Did not explicity match on previous vote. Checking for balance post-matching. 

pres <- read_csv("pres_county_80-16.csv") # Presidential Pre-trend data 
pres$FIPS <- ((pres$FIPS1*1000)+ pres$FIPS2) # create similar FIPS code 
matchedhospitalcounties <- merge(matchedhospitalcounties, pres, by="FIPS")
matchedhospitalcountiesr <- merge(matchedhospitalcountiesr, pres, by="FIPS")
matchedhospitalcountiesrns <- merge(matchedhospitalcountiesrns, pres, by="FIPS")
matchedhospitalcountiesns <- merge(matchedhospitalcountiesns, pres, by="FIPS")
# Pre-treatment Partisan Balance (Primary Match)
summary(lm(matchedhospitalcounties$dm2000~  matchedhospitalcounties$treatedever))
summary(lm(matchedhospitalcounties$dm2004~ matchedhospitalcounties$treatedever))

treatedc <- matchedhospitalcounties[which(matchedhospitalcounties$treatedever==1),]
untreatedc <- matchedhospitalcounties[which(matchedhospitalcounties$treatedever==0),]

weighted.mean(treatedc$dm2000, treatedc$Population_2010)
weighted.mean(untreatedc$dm2000, untreatedc$Population_2010)

weighted.mean(treatedc$dm2004, treatedc$Population_2010)
weighted.mean(untreatedc$dm2004, untreatedc$Population_2010)

# Pre-treatment Partisan Balance (Primary Match + Texas Robustness)
summary(lm(matchedhospitalcountiesr$dm2000~ matchedhospitalcountiesr$treatedever2))
summary(lm(matchedhospitalcountiesr$dm2004~ matchedhospitalcountiesr$treatedever2))

treatedcr <- matchedhospitalcountiesr[which(matchedhospitalcountiesr$treatedever2==1),]
untreatedcr <- matchedhospitalcountiesr[which(matchedhospitalcountiesr$treatedever2==0),]

weighted.mean(treatedcr$dm2000, treatedcr$Population_2010)
weighted.mean(untreatedcr$dm2000, untreatedcr$Population_2010)

weighted.mean(treatedcr$dm2004, treatedcr$Population_2010)
weighted.mean(untreatedcr$dm2004, untreatedcr$Population_2010)



# Pre-treatment Partisan Balance (In State Match)
summary(lm(matchedhospitalcountiesns$dm2000~ matchedhospitalcountiesns$treatedever))
summary(lm(matchedhospitalcountiesns$dm2004~ matchedhospitalcountiesns$treatedever))


treatedcns <- matchedhospitalcountiesns[which(matchedhospitalcountiesns$treatedever2==1),]
untreatedns <- matchedhospitalcountiesns[which(matchedhospitalcountiesns$treatedever2==0),]

weighted.mean(treatedcns$dm2000, treatedcns$Population_2010)
weighted.mean(untreatedns$dm2000, untreatedns$Population_2010)

weighted.mean(treatedcns$dm2004, treatedcns$Population_2010)
weighted.mean(untreatedns$dm2004, untreatedns$Population_2010)


# Pre-treatment Partisan Balance (In State Match + Texas Robustness)
summary(lm(matchedhospitalcountiesrns$dm2000~ matchedhospitalcountiesrns$treatedever2))
summary(lm(matchedhospitalcountiesrns$dm2004~ matchedhospitalcountiesrns$treatedever2))

# No Statistically Significant Differences. 

######################################################################
#               Loading and Cleaning CES Data
#####################################################################
ces <- read_dta("cumulative_2006-2022.dta")
cespolicy <- read_dta("cumulative_ces_policy_preferences.dta")

ces$female[ces$gender==1]<- 0
ces$female[ces$gender==2]<- 1

ces$nonwhite[ces$race==1] <- 0
ces$nonwhite[ces$race!=1] <- 1

ces$republican <- NA
ces$republican[ces$pid7==5] <- 1 
ces$republican[ces$pid7==6] <- 2 
ces$republican[ces$pid7==7] <- 3
ces$republican[ces$pid7==4] <- 0
ces$republican[ces$pid7==3] <- -1
ces$republican[ces$pid7==2] <- -2
ces$republican[ces$pid7==1] <- -3


ces$faminc[ces$faminc==13] <- NA 
ces$faminc[ces$faminc==14] <- NA 

ces$unemployed <- 0 
ces$unemployed[ces$employ==4] <- 1
ces$unemployed[ces$employ==3] <- 1

ces$christian <- NA 
ces$christian[ces$religion>=5] <- 0 
ces$christian[ces$religion<5] <- 1 

ces$govapproval <- NA 
ces$govapproval[ces$approval_gov==1] <- 3
ces$govapproval[ces$approval_gov==2] <- 2
ces$govapproval[ces$approval_gov==3] <- 1
ces$govapproval[ces$approval_gov==4] <- 0

ces$presapproval <- NA 
ces$presapproval[ces$approval_pres==1] <- 3
ces$presapproval[ces$approval_pres==2] <- 2
ces$presapproval[ces$approval_pres==3] <- 1
ces$presapproval[ces$approval_pres==4] <- 0


ces$voted_gop < - NA 
ces$voted_gop[ces$voted_pres_party==2] <- 1
ces$voted_gop[ces$voted_pres_party==1] <- 0


ces$gopgov <- NA 
ces$gopgov[ces$voted_gov==1] <- 0
ces$gopgov[ces$voted_gov==2] <- 1


# Merge with Geographic Data 
ces_full <- merge(ces, rural, by.x = "county_fips", by.y ="FIPS")

ces_full$nonexpansion <- 0 
ces_full$nonexpansion[ces_full$State=="TX" |
                        ces_full$State=="TN" |
                        ces_full$State=="FL" |
                        ces_full$State=="NC" |
                        ces_full$State=="GA" |  ces_full$State=="AL" | ces_full$State=="MS" | 
                        ces_full$State=="KS" | ces_full$State=="WY" | ces_full$State=="WI" | 
                        ces_full$State=="SC" | ces_full$State=="SD" | ces_full$State=="OK" | ces_full$State=="MO" | ces_full$State=="NE" | ces_full$State=="ID" |
                        ces_full$State=="ME"] <- 1 



# Merge with policy ID 
policy_full  <- inner_join(ces_full, cespolicy)
# Subset to White Respondents 
policy_full <- policy_full[which(policy_full$nonwhite!=1),]
# All Respondents 
policy_full2  <- inner_join(ces_full, cespolicy)


# ACA Support (White Respondents)
policy_full$aca_support <- NA 
policy_full$aca_support[policy_full$healthcare_aca==1] <- 0
policy_full$aca_support[policy_full$healthcare_aca==2] <- 1

# 2014 Needs to be flipped (due to apparent CES coding error in version)
policy_full$aca_support[policy_full$healthcare_aca==1 & policy_full$year==2014] <- 1
policy_full$aca_support[policy_full$healthcare_aca==2 & policy_full$year==2014] <- 0



# ACA Support (All Respondents)
policy_full2$aca_support <- NA 
policy_full2$aca_support[policy_full2$healthcare_aca==1] <- 0
policy_full2$aca_support[policy_full2$healthcare_aca==2] <- 1

# 2014 Needs to be flipped (due to apparent CES coding error in version)
policy_full2$aca_support[policy_full2$healthcare_aca==1 & policy_full2$year==2014] <- 1
policy_full2$aca_support[policy_full2$healthcare_aca==2 & policy_full2$year==2014] <- 0


# Demonstrates that "flip" was necessary. 

ruralplot <- policy_full %>%
  group_by(rural, year) %>%
  summarize(meanaca = weighted.mean(aca_support, w= weight_cumulative,  na.rm=T))

ruralplot <- ruralplot[complete.cases(ruralplot$meanaca),]

ruralplot$group <- NA 
ruralplot$group[ruralplot$rural==1] <- "Rural"
ruralplot$group[ruralplot$rural==0] <- "Not Rural"


ruralplot$meanaca <- ruralplot$meanaca*100

c <- ggplot(ruralplot, aes(x=year, y= meanaca, color=group, fill=group)) +   
  geom_line(size=1)   + ylim(20,80) + 
  labs(y="ACA Support %", 
       x = "Year", 
       color="Community", 
       fill="Community") + scale_color_manual(values=c("gray","black")) + scale_fill_manual(values=c("gray","black")) + 
  theme_classic()


######################################################################
#              Merging Matched County &  CES Data
#####################################################################

# Merge with Primary Matched County data 
ces_hosp <- merge(ces_full, m.data, by.x = "county_fips", by.y = "FIPS")
# Merge with Primary Matched County + Texas Robust data 
ces_hospr <- merge(ces_full, m.datar, by.x = "county_fips", by.y = "FIPS")
# Merge with In-State Matched data 
ces_hosp_instate <- merge(ces_full, m.out0_statematch, by.x = "county_fips", by.y = "FIPS")
# Merge with In-State + Texas Robust Data
ces_hosp_instater <- merge(ces_full, m.datars, by.x = "county_fips", by.y = "FIPS")


# White Respondents 
ces_hosp2 <- merge(ces_hosp, policy_full, by="case_id")
ces_hosp2r <- merge(ces_hospr, policy_full, by="case_id")
ces_hosp_instate2 <- merge(ces_hosp_instate, policy_full, by="case_id")
ces_hosp_instate2r <- merge(ces_hosp_instater, policy_full, by="case_id")

# All Respondents 
ces_hosp2a <- merge(ces_hosp, policy_full2, by="case_id")
ces_hosp2ra <- merge(ces_hospr, policy_full2, by="case_id")
ces_hosp_instate2a <- merge(ces_hosp_instate, policy_full2, by="case_id")
ces_hosp_instate2ra <- merge(ces_hosp_instater, policy_full, by="case_id")

### Creating Treatment Variables 

# White Respondents 
# Treated Primary Matched County data 
ces_hosp2$treated <- 0
ces_hosp2$treated[ces_hosp2$year.x >= ces_hosp2$`Closure Year` &  ces_hosp2$treatedever==1] <- 1 
# Treated In State Match 
ces_hosp_instate2$treated <- 0
ces_hosp_instate2$treated[ces_hosp_instate2$year.x >= ces_hosp_instate2$`Closure Year` &  ces_hosp_instate2$treatedever==1] <- 1 
## Treated Primary Matched County data  + Texas Case 
ces_hosp2r$treated <- 0
ces_hosp2r$treated[ces_hosp2r$year.x>=ces_hosp2r$`Closure Year`  &  ces_hosp2r$treatedever2==1] <- 1 
# Treated In State Match + Texas Case 
ces_hosp_instate2r$treated  <- 0 
ces_hosp_instate2r$treated[ces_hosp_instate2r$year.x>=ces_hosp_instate2r$`Closure Year` & ces_hosp_instate2r$treatedever2==1] <- 1 

# All Respondents 
# Treated Primary Matched County data 
ces_hosp2a$treated <- 0
ces_hosp2a$treated[ces_hosp2a$year.x >= ces_hosp2a$`Closure Year` &  ces_hosp2a$treatedever==1] <- 1 
# Treated In State Match 
ces_hosp_instate2a$treated <- 0
ces_hosp_instate2a$treated[ces_hosp_instate2a$year.x >= ces_hosp_instate2a$`Closure Year` &  ces_hosp_instate2a$treatedever==1] <- 1 
## Treated Primary Matched County data  + Texas Case 
ces_hosp2ra$treated <- 0
ces_hosp2ra$treated[ces_hosp2ra$year.x>=ces_hosp2ra$`Closure Year`  &  ces_hosp2ra$treatedever2==1] <- 1 
# Treated In State Match + Texas Case 
ces_hosp_instate2ra$treated  <- 0 
ces_hosp_instate2ra$treated[ces_hosp_instate2ra$year.x>=ces_hosp_instate2ra$`Closure Year` & ces_hosp_instate2ra$treatedever2==1] <- 1 



######################################################################
#            Probing Individual-Level Balance
#####################################################################

# Without Milam County 
descripstreated <- ces_hosp2a[which(ces_hosp2a$treatedever==1),]
descripsuntreated <- ces_hosp2a[which(ces_hosp2a$treatedever==0),]

weighted.mean(ces_hosp2a$age.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$age.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$age.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$republican.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$republican.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$republican.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$female.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$female.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$female.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$educ.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$educ.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$educ.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$faminc.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$faminc.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$faminc.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$christian.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$christian.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$christian.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$gopgov.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(ces_hosp2a$voted_gop.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)



# With Milam County
descripstreated <- ces_hosp2ra[which(ces_hosp2ra$treatedever2==1),]
descripsuntreated <- ces_hosp2ra[which(ces_hosp2ra$treatedever2==0),]

weighted.mean(ces_hosp2a$age.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$age.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$age.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$republican.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$republican.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$republican.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$female.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$female.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$female.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$educ.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$educ.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$educ.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$faminc.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$faminc.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$faminc.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$christian.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(descripstreated$christian.x, w=descripstreated$weight_cumulative.x, na.rm = T)
weighted.mean(descripsuntreated$christian.x, w=descripsuntreated$weight_cumulative.x, na.rm = T)

weighted.mean(ces_hosp2a$gopgov.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)
weighted.mean(ces_hosp2a$voted_gop.x, w=ces_hosp2a$weight_cumulative.x, na.rm = T)


######################################################################
#            Regression Analyses (Nationwide Matches)
#####################################################################
ces_hosp2_nofips <- ces_hosp2[which(ces_hosp2$county_fips.x!="NA"),]   # drop obvs with no FIPS Code due to error clusters if any 

#### Matching Main Results without Milam (White Respondents)
matching1 <- lm(voted_gop.x~  treated +  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2 <- lm(gopgov.x~  treated + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching1c <- lm(voted_gop.x~  treated+ republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2c <- lm(gopgov.x~  treated +republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

matching1ci <- lm(voted_gop.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2ci <- lm(gopgov.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

match1cla <- coeftest(matching1 , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2cla <- coeftest(matching2, vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match1clb <- coeftest(matching1c , vcov = vcovCL, cluster = ~county_fips.x+State.x)
match2clb <- coeftest(matching2c, vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match1clc <- coeftest(matching1ci , vcov = vcovCL, cluster = ~county_fips.x+State.x)
match2clc <- coeftest(matching2ci, vcov = vcovCL, cluster = ~county_fips.x+State.x) 


# main text matching results 
stargazer(matching1, matching1c, matching1ci) 
stargazer(match1cla, match1clb, match1clc) 
stargazer(matching2, matching2c, matching2c) 
stargazer(match2cla, match2clb, match2clc) 



#### Within Matched Pair Presidential Regression 
base1sc <- lm(voted_gop.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)
#### Within County Presidential Regression 
base1fc <- lm(voted_gop.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 
base2c_cl <- coeftest(base1fc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1sc, base1fc)
stargazer(base1c_cl, base2c_cl)


#### Within Matched Pair/County Gov Results
base1c <- lm(gopgov.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)
base1sc <- lm(gopgov.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)


base1c_cl <- coeftest(base1c, vcov = vcovCL, cluster = ~county_fips.x+State) 
base1sc_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1c, base1sc )
stargazer(base1c_cl, base1sc_cl)



#### Matching Results without Milam (All Respondents)
ces_hosp2_nofips <- ces_hosp2a[which(ces_hosp2a$county_fips.x!="NA"),] # drop obvs with no FIPS Code due to error clusters 

matching1 <- lm(voted_gop.x~  treated +  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2 <- lm(gopgov.x~  treated + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching1c <- lm(voted_gop.x~  treated+ republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2c <- lm(gopgov.x~  treated +republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

matching1ci <- lm(voted_gop.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2ci <- lm(gopgov.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

match1cla <- coeftest(matching1 , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2cla <- coeftest(matching2, vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match1clb <- coeftest(matching1c , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2clb <- coeftest(matching2c, vcov = vcovCL, cluster = ~county_fips.x+State.x)  
match1clc <- coeftest(matching1ci , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2clc <- coeftest(matching2ci, vcov = vcovCL, cluster = ~county_fips.x+State.x)


# main text matching results 
stargazer(matching1, matching1c, matching1ci) 
stargazer(match1cla, match1clb, match1clc) 
stargazer(matching2, matching2c, matching2c) 
stargazer(match2cla, match2clb, match2clc) 


#### Within Matched Pair Presidential Regression 
base1sc <- lm(voted_gop.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)
####  Within County Presidential Regression
base1fc <- lm(voted_gop.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 
base2c_cl <- coeftest(base1fc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1sc, base1fc)
stargazer(base1c_cl,  base2c_cl)


#### Within Matched Pair/County Gov Results
base1c <- lm(gopgov.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)
base1sc <- lm(gopgov.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1c, vcov = vcovCL, cluster = ~county_fips.x+State) 
base1sc_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1c, base1sc)
stargazer(base1c_cl, base1sc_cl)



#### Matching Main Results With Milam Texas Case (White Respondents)
ces_hosp2_nofips <- ces_hosp2r[which(ces_hosp2r$county_fips.x!="NA"),] # drop obvs with no FIPS Code due to error clusters 

matching1 <- lm(voted_gop.x~  treated +  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2 <- lm(gopgov.x~  treated + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching1c <- lm(voted_gop.x~  treated+ republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2c <- lm(gopgov.x~  treated +republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

matching1ci <- lm(voted_gop.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2ci <- lm(gopgov.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

match1cla <- coeftest(matching1 , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2cla <- coeftest(matching2, vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match1clb <- coeftest(matching1c , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2clb <- coeftest(matching2c, vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match1clc <- coeftest(matching1ci , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2clc <- coeftest(matching2ci, vcov = vcovCL, cluster = ~county_fips.x+State.x) 


# main text matching results 
stargazer(matching1, matching1c, matching1ci) 
stargazer(match1cla, match1clb, match1clc) 
stargazer(matching2, matching2c, matching2c) 
stargazer(match2cla, match2clb, match2clc) 


#### Within Matched Pair Presidential Regression 
base1sc <- lm(voted_gop.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)
####  Within County Presidential Regression
base1fc <- lm(voted_gop.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 
base2c_cl <- coeftest(base1fc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1sc, base1fc)
stargazer(base1c_cl,  base2c_cl)


#### Within Matched Pair/County Gov Results
base1c <- lm(gopgov.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)
base1sc <- lm(gopgov.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1c, vcov = vcovCL, cluster = ~county_fips.x+State) 
base1sc_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1c, base1sc)
stargazer(base1c_cl, base1sc_cl)



#### Matching Main Results With Milam Texas Case (All Respondents)
ces_hosp2_nofips <- ces_hosp2ra[which(ces_hosp2ra$county_fips.x!="NA"),] # drop obvs with no FIPS Code due to error clusters 

matching1 <- lm(voted_gop.x~  treated +  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2 <- lm(gopgov.x~  treated + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching1c <- lm(voted_gop.x~  treated+ republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2c <- lm(gopgov.x~  treated +republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                   faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

matching1ci <- lm(voted_gop.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x+  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2ci <- lm(gopgov.x~  treated*republican.x + age.x + female.x + educ.x + nonexpansion + christian.x + 
                    faminc.x + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

match1cla <- coeftest(matching1 , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2cla <- coeftest(matching2, vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match1clb <- coeftest(matching1c , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2clb <- coeftest(matching2c, vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match1clc <- coeftest(matching1ci , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2clc <- coeftest(matching2ci, vcov = vcovCL, cluster = ~county_fips.x+State.x)


# main text matching results 
stargazer(matching1, matching1c, matching1ci) 
stargazer(match1cla, match1clb, match1clc) 
stargazer(matching2, matching2c, matching2c) 
stargazer(match2cla, match2clb, match2clc) 


#### Within Matched Pair Presidential Regression 
base1sc <- lm(voted_gop.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)
####  Within County Presidential Regression
base1fc <- lm(voted_gop.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 
base2c_cl <- coeftest(base1fc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1sc, base1fc)
stargazer(base1c_cl,  base2c_cl)


#### Within Matched Pair/County Gov Results
base1c <- lm(gopgov.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)
base1sc <- lm(gopgov.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1c, vcov = vcovCL, cluster = ~county_fips.x+State) 
base1sc_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1c, base1sc)
stargazer(base1c_cl, base1sc_cl)


##### MECHANISMS 
ces_hosp2_p <- merge(ces_hosp2_nofips, policy_full, by="case_id") # White Respondents 
ces_hosp2_pa <- merge(ces_hosp2_nofips, policy_full2, by="case_id") # All Respondents 

obamayear <- ces_hosp2_p[which(ces_hosp2_p $year.y<=2016),]
trumpyear<- ces_hosp2_p[which(ces_hosp2_p $year.y>=2017),]


obamayeara <- ces_hosp2_pa[which(ces_hosp2_pa$year.y<=2016),]
trumpyeara<- ces_hosp2_pa[which(ces_hosp2_pa$year.y>=2017),]


presapp1 <- lm(presapproval  ~  treated + distance  + age.x + female.x + educ.x + christian.x + 
                 faminc.x, data=obamayear, weights = weight_cumulative.x)
presapp2 <- lm(presapproval  ~  treated + distance  + age.x + female.x + educ.x + christian.x + 
                 faminc.x, data=trumpyear, weights = weight_cumulative.x)

aca1 <- lm(aca_support.y ~  treated  + distance + age.x + female.x + educ.x + christian.x + 
             faminc.x, data=obamayear, weights = weight_cumulative.x)
aca2 <- lm(aca_support.y ~  treated  + distance  + age.x + female.x + educ.x  + christian.x + 
             faminc.x, data=trumpyear, weights = weight_cumulative.x)

govsapp1 <- lm(govapproval ~  treated  + distance + age.x + female.x + educ.x  + christian.x + 
                 faminc.x, data=obamayear, weights = weight_cumulative.x)
govsapp2 <- lm(govapproval ~  treated + distance  + age.x + female.x + educ.x  + christian.x + 
                 faminc.x, data=trumpyear, weights = weight_cumulative.x)

stargazer(presapp1, presapp2, aca1, aca2, govsapp1, govsapp2)


presapp1 <- lm(presapproval  ~  treated  + distance + age.x + female.x + educ.x  + christian.x + 
                 faminc.x + nonwhite, data=obamayeara, weights = weight_cumulative.x)
presapp2 <- lm(presapproval  ~  treated + distance  + age.x + female.x + educ.x  + christian.x + 
                 faminc.x + nonwhite, data=trumpyeara, weights = weight_cumulative.x)

aca1 <- lm(aca_support.y ~  treated + distance + age.x + female.x + educ.x  + christian.x + 
             faminc.x + nonwhite, data=obamayeara, weights = weight_cumulative.x)
aca2 <- lm(aca_support.y ~  treated + distance + age.x + female.x + educ.x  + christian.x + 
             faminc.x + nonwhite, data=trumpyeara, weights = weight_cumulative.x)

govsapp1 <- lm(govapproval ~  treated + distance + age.x + female.x + educ.x + christian.x + 
                 faminc.x + nonwhite, data=obamayeara, weights = weight_cumulative.x)
govsapp2 <- lm(govapproval ~  treated + distance  + age.x + female.x + educ.x  + christian.x + 
                 faminc.x + nonwhite, data=trumpyeara, weights = weight_cumulative.x)

stargazer(presapp1, presapp2, aca1, aca2, govsapp1, govsapp2)


######################################################################
#   Regression Analyses (Within State Matches Robustness )
#####################################################################
#### Matching Within State Results (All Respondents)
ces_hosp2_nofips <- ces_hosp_instate2a[which(ces_hosp_instate2a$county_fips.x!="NA"),]   # drop obvs with no FIPS Code due to error clusters if any 

matching1 <- lm(voted_gop.x~  treated +  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2 <- lm(gopgov.x~  treated + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

match1cla <- coeftest(matching1 , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2cla <- coeftest(matching2, vcov = vcovCL, cluster = ~county_fips.x+State.x) 

# Within State matching results 
stargazer(matching1, matching2) 
stargazer(match1cla, match2cla) 



#### Within Matched Pair  Within State Presidential Regression 
base1sc <- lm(voted_gop.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)
#### Within County Within State Presidential Regression 
base1fc <- lm(voted_gop.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 
base2c_cl <- coeftest(base1fc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1sc, base1fc)
stargazer(base1c_cl, base2c_cl)


#### Within Matched Pair/County  Within State Gov Results
base1c <- lm(gopgov.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)
base1sc <- lm(gopgov.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)


base1c_cl <- coeftest(base1c, vcov = vcovCL, cluster = ~county_fips.x+State) 
base1sc_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1c, base1sc )
stargazer(base1c_cl, base1sc_cl)


# Treated In State Match + Texas Case All Respondents 
ces_hosp2_nofips <- ces_hosp_instate2ra[which(ces_hosp_instate2ra$county_fips.x!="NA"),]   # drop obvs with no FIPS Code due to error clusters if any 

matching1 <- lm(voted_gop.x~  treated +  distance + as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)
matching2 <- lm(gopgov.x~  treated + distance +  as.factor(year.y), data=ces_hosp2_nofips, weights = weight_cumulative.x)

match1cla <- coeftest(matching1 , vcov = vcovCL, cluster = ~county_fips.x+State.x) 
match2cla <- coeftest(matching2, vcov = vcovCL, cluster = ~county_fips.x+State.x) 


# Within State matching results 
stargazer(matching1, matching2) 
stargazer(match1cla, match2cla) 

#### Within Matched Pair  Within State Presidential Regression 
base1sc <- lm(voted_gop.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)
#### Within County Within State Presidential Regression 
base1fc <- lm(voted_gop.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)

base1c_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 
base2c_cl <- coeftest(base1fc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1sc, base1fc)
stargazer(base1c_cl, base2c_cl)


#### Within Matched Pair/County  Within State Gov Results
base1c <- lm(gopgov.x ~  treated +  as.factor(county_fips.x), data=ces_hosp2_nofips, weights = weight_cumulative.x)
base1sc <- lm(gopgov.x ~  treated +  as.factor(subclass), data=ces_hosp2_nofips, weights = weight_cumulative.x)


base1c_cl <- coeftest(base1c, vcov = vcovCL, cluster = ~county_fips.x+State) 
base1sc_cl <- coeftest(base1sc, vcov = vcovCL, cluster = ~county_fips.x+State) 

stargazer(base1c, base1sc )
stargazer(base1c_cl, base1sc_cl)




######################################################################
#               Placebo Robustness Analyses 
#####################################################################
ces_hosp2$treatwindow <- NA 
ces_hosp2$treatwindow <- ces_hosp2$year.x - ces_hosp2$`Closure Year`

# One year 
ces_hosp2$placebo <- NA
ces_hosp2$placebo[ces_hosp2$treatwindow<=-1] <- 1 
ces_hosp2$placebo[ces_hosp2$treatwindow>=0] <- 0 
ces_hosp2$placebo[ces_hosp2$treatedever==0] <- 0 
# two years 
ces_hosp2$placebo2a <-NA
ces_hosp2$placebo2a[ces_hosp2$treatwindow<=-2] <- 1
ces_hosp2$placebo2a[ces_hosp2$treatwindow>=0] <- 0 
ces_hosp2$placebo2a[ces_hosp2$treatedever==0] <- 0 
# three years 
ces_hosp2$placebo3 <-NA
ces_hosp2$placebo3[ces_hosp2$treatwindow<=-3] <- 1
ces_hosp2$placebo3[ces_hosp2$treatwindow>=0] <- 0 
ces_hosp2$placebo3[ces_hosp2$treatedever==0] <- 0 

# Four Years
ces_hosp2$placebo2 <-NA
ces_hosp2$placebo2[ces_hosp2$treatwindow<=-4] <- 1
ces_hosp2$placebo2[ces_hosp2$treatwindow>=0] <- 0 
ces_hosp2$placebo2[ces_hosp2$treatedever==0] <- 0 


p1 <- lm(voted_gop.x~  placebo + republican.x+ age.x + female.x + educ.x + nonexpansion + christian.x + 
           faminc.x + `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
           `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + as.factor(State)+ as.factor(year.y), data=ces_hosp2, weights = weight_cumulative.x)


p2 <- lm(voted_gop.x~  placebo2a + republican.x+ age.x + female.x + educ.x + nonexpansion + christian.x + 
           faminc.x + `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
           `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + as.factor(State)+ as.factor(year.y), data=ces_hosp2, weights = weight_cumulative.x)

p3 <- lm(voted_gop.x~  placebo3 + republican.x+ age.x + female.x + educ.x + nonexpansion + christian.x + 
           faminc.x + `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
           `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + as.factor(State)+ as.factor(year.y), data=ces_hosp2, weights = weight_cumulative.x)

p4<- lm(voted_gop.x~  placebo2 + republican.x+ age.x + female.x + educ.x + nonexpansion + christian.x + 
          faminc.x + `% Rural raw value` + `Uninsured adults raw value`+ `% 65 and older raw value` +`% Non-Hispanic African American raw value` + 
          `% Hispanic raw value` +`Median household income raw value`  + `Poor or fair health raw value` + as.factor(State)+ as.factor(year.y), data=ces_hosp2, weights = weight_cumulative.x)


stargazer(p1, p2, p3, p4, 
          dep.var.labels=c("Voted Republican President"),
          covariate.labels=c("Placebo (t-1)","Placebo (t-2)", "Placebo (t-3)", "Placebo (t-4)",
                             "Partisanship",
                             "Age", "Female", "Education", "Non-expansion State", "Christian", "Income", 
                             "Rural", "Uninsured %", "Over 65%", "% Black", "% Hispanic", "Median Income" , "% Poor or Fair Health"), omit=c("State", "year"), 
          style = "ajps", 
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 





######################################################################
#         County Election Data Robustness Analysis 
#####################################################################
pres2<- read_csv("countypres_2000-2020.csv") # County Election Data 2000-2020 
# Cleaning and Rearranging Data for Merge 
pres2dem <- pres2[which(pres2$party=="DEMOCRAT"),]
pres2gop <- pres2[which(pres2$party=="REPUBLICAN"),]

pres2dem$dt<- pres2dem$candidatevotes
pres2gop$rt <- pres2gop$candidatevotes

keep <- c("year", "county_fips", "dt")

pres2dem <- pres2dem[keep]

prescountyfull <- merge(pres2dem, pres2gop, by=c("county_fips", "year"))

prescountyfull$goptwoparty <- ((prescountyfull$rt)/(prescountyfull$rt+prescountyfull$dt))
prescountyfull$gopoverall <- ((prescountyfull$rt)/(prescountyfull$rt+prescountyfull$totalvotes))

#With matched counties (No Milam) 
prescountymatched <- merge(prescountyfull, matchedhospitalcounties, by.x = "county_fips", by.y = "FIPS")

prescountymatched$treated <- 0 
prescountymatched$treated[prescountymatched$`Closure Year`<=prescountymatched$year] <- 1 


pooled1 <- lm(prescountymatched$goptwoparty ~ prescountymatched$treated + as.factor(prescountymatched$county_fips))
pooled2 <- lm(prescountymatched$goptwoparty ~ prescountymatched$treated + prescountymatched$`% 65 and older raw value` +
                prescountymatched$`Median household income raw value` + prescountymatched$`% Non-Hispanic African American raw value` +
                prescountymatched$`% Non-Hispanic African American raw value` + prescountymatched$`Uninsured adults raw value` + as.factor(prescountymatched$state))


pooled1cl <- coeftest(pooled1, vcov = vcovCL, cluster = ~prescountymatched$county_fips + prescountymatched$state) # one tail significant, close to single
pooled2cl <- coeftest(pooled2, vcov = vcovCL, cluster = ~prescountymatched$county_fips + prescountymatched$state) # one tail significant, close to single


stargazer(pooled1,  pooled2)
stargazer(pooled1cl,pooled2cl)



#With matched counties (With Milam) 
prescountymatchedr <- merge(prescountyfull, matchedhospitalcountiesr, by.x = "county_fips", by.y = "FIPS")

prescountymatchedr$treated <- 0 
prescountymatchedr$treated[prescountymatchedr$`Closure Year`<=prescountymatchedr$year] <- 1 


pooled1 <- lm(prescountymatchedr$goptwoparty ~ prescountymatchedr$treated + as.factor(prescountymatchedr$county_fips))
pooled2 <- lm(prescountymatchedr$goptwoparty ~ prescountymatchedr$treated + prescountymatchedr$`% 65 and older raw value` +
                prescountymatchedr$`Median household income raw value` + prescountymatchedr$`% Non-Hispanic African American raw value` +
                prescountymatchedr$`% Non-Hispanic African American raw value` + prescountymatchedr$`Uninsured adults raw value` + as.factor(prescountymatchedr$state))


pooled1cl <- coeftest(pooled1, vcov = vcovCL, cluster = ~prescountymatchedr$county_fips + prescountymatchedr$state) # one tail significant, close to single
pooled2cl <- coeftest(pooled2, vcov = vcovCL, cluster = ~prescountymatchedr$county_fips + prescountymatchedr$state) # one tail significant, close to single


stargazer(pooled1,  pooled2)
stargazer(pooled1cl,pooled2cl)



